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Abstract. We calculate the stability of a multiple relaxation region MHD 
(MRXMHD) plasma, or stepped-Beltrami plasma, using both variational and tearing 
mode treatments. The configuration studied is a periodic cylinder. In the 
variational treatment, the problem reduces to an eigenvalue problem for the interface 
displacements. For the tearing mode treatment, analytic expressions for the tearing 
mode stability parameter A', being the jump in the logarithm in the helical flux across 
the resonant surface, are found. The stability of these treatments is compared for m = 1 
displacements of an illustrative RFP-like conflguration, comprising two distinct plasma 
regions. For pressure- less conflgurations, we flnd the marginal stability conclusions 
of each treatment to be identical, conflrming analytic results in the literature. The 
tearing mode treatment also resolves ideal MHD unstable solutions for which A' ^ oo: 
these correspond to displacement of a resonant interface. Wall stabilisation scans 
resolve the internal and external ideal kink. Scans with increasing pressure are also 
performed: these indicate that both variational and tearing mode treatments have the 
same stability trends with and show pressure stabilisation in conflgurations with 
increasing edge pressure. Combined, our results suggest that MRXMHD conflgurations 
which are stable to ideal perturbations plus tearing modes are automatically in a stable 
state. Such configurations, and their stability properties, are of emerging importance 
in the quest to find mathematically rigorous solutions of ideal MHD force balance in 
3D geometry. 
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1. Introduction 

Recently, Hole et al [1] proposed a model for a partially relaxed plasma-vacuum system. 
The purpose of the model, which abandons all but a small number of flux surfaces, 
is to provide a mathematically rigorous foundation for ideal MHD equilibria in 3D 
configurations. The model appeals to both chaotic field lines, that flatten the pressure 
gradient in chaotic regions, and Taylor relaxation, which force the plasma gradient to be 
zero in Taylor-relaxed regions. The model consists of a stepped pressure profile, where 
the steps correspond to ideal MHD barriers across which can be supported a pressure 
or field jump, or a jump in rotational transform. Our overarching objective is the 
development of an equilibrium solver for 3D plasmas built on a stepped pressure profile 
model. In the 3D case, we envisage that the barriers can be chosen to be non-resonant 
KAM surfaces that survive the onset of field line chaos intrinsic to 3D equilibria. In 
between the interfaces, the field is Beltrami, such that V x B = yuB. The boundary 
condition across the interfaces is the continuity of total pressure p + B^/2fio. Such a 
model, which we term a MRXMHD (multiple relaxation regions MHD) model, raises a 
number of questions. How should the equilibrium be constrained? How much jump in 
pressure and/or rotational transform * can each interface support? Are the interfaces 
stable to deformation? Can the class of stability shed information onto other quasi- 
relaxed phenomena? 

Previous work has focused on the equilibrium constraints [2, 3], construction 
of a numerical algorithm for calculation of Beltrami fields between interfaces in 3D 
configurations [4], and a variational principle for the equilibrium and stability of the 
multiple interface configuration in cylindrical plasmas [1]. We have also explored the 
relationship between relaxed plasma equilibrium models discussed here, and entropy 
related plasma self-organisation principles [5]. The motivation of this paper is to 
understand the nature of MRXMHD modes identified from the variational principle. 
Our work complements a separate in-press publication [6] that unifies relaxed and ideal 
MHD principles for constructing global solutions comprising mixed relaxed and ideal 
regions. 

Recently, Tassi et al [7], performed a tearing mode stability treatment on stepped 
/i force-free equilibria close to Taylor relaxed states. Their motivation was to develop 
a mechanism for the formation of cyclic Quasi-Single-Helicity (QSH) states observed in 
Reverse Field Pinches [8]. They consider a cylindrical plasma divided into two different 
Beltrami regions, and encased in a perfectly conducting shell, and compute the tearing 
mode stability parameter A' at a resonant radius for a helical fiux perturbation 
Xi(r) = mBziir) — KrBgi{r). Here, r is the radial coordinate, and m and k the poloidal 
and axial wave number. Tassi et al find critical values of the jump in /x, beyond which 
the RFP-like plasma is unstable. Based on these, they postulate the QSH state may be 
viewed as a small, cyclic departure from a Taylor-relaxed state. 

In this work, we extend the tearing mode stability treatment of Tassi et al [7] to 
plasmas with finite pressure and a vacuum region, and compare stability conclusions 



of our variational treatment to that of a tearing mode stability analysis. Our paper is 
arranged as follows : Sec. 2 summarises the variational model of stepped pressure 
profile plasmas, presented in Hole et al [1], and introduces a tearing mode model. 
Section 3 treats MRXMHD plasmas in cylindrical geometry, yielding stability parameter 
expressions for both the variational and tearing mode treatments. In Sec. 4, we 
compute stability for an example configuration, draw comparisons between the stability 
conclusions based on variational and tearing mode treatments, and explore marginal 
stability limits in wave-number space as a function of pressure. Finally, Sec. 5 contains 
concluding remarks. 

2. Multiple-interface plasma-vacuum model 

The system comprises Taylor-relaxed plasma regions, each separated by an ideal 
MHD barrier. The outermost plasma region is enclosed by a vacuum, and encased in a 
perfectly conducting wall. Figure 1(a) shows the geometry of the system, and introduces 
the nomenclature used to describe the region and interfaces. The regions TZi comprise 
the plasma regions IZi = Vi, ...,71^ = Vn and the vacuum region T^at+i = V. Each 
plasma region Vi is bounded by the inner and outer ideal MHD interfaces Xj_i, and Xj 
respectively, whilst the vacuum is encased by the perfectly conducting wall W. 




Figure 1. Schematic of magnetic geometry (a), showing ideal MHD barriers 
Xj, the conducting wall W, plasma regions Vi and the vacuum V. Panel (b) 
shows the q profile used for stability studies in Sec. 4, with //i = 2 (core) and 
//2 = 3.6 (edge). 

2.1. A variational description 

In previous work [1] we outlined our variational principle, which lies between that of 
Kruskal &; Kulsrud [9] — minimization of total energy W = J[B^ /2 + p/ (7 — 1)] (where 
p is plasma pressure and 7 the ratio of specific heats) under the uncountable infinity of 



constraints provided by applying ideal MHD within each fluid element — and the relaxed 
MHD of Woltjer [10] and Taylor [11] — minimization of W holding only the two global 
toroidal and poloidal magnetic fluxes, and the single global ideal-MHD helicity invariant 
H = J A ■ B, constant. In summary, the energy functional could be written 

N N N 

i=l 1=1 1=1 

where /ij and Ui are Lagrange multipliers, and 

U^= [ dr'(^ + ^], (2) 



7^,: 



7-1 2/io 



M,= / rfr^p//^, (3) 

Hi = [ dr^A ■ V X A + / dl ■ A / dl ■ A - / dl ■ A / dl ■ A. (4) 

The term f/j is the potential energy, Mj the plasma mass, and Hi the magnetic helicity 
in each region TZi. In Eqs. (2) - (4), dr^ is a volume element, 7 the ratio of speciflc 
heats, and Pi, Bi and Aj the equilibrium pressure, magnetic held strength and vector 
potential respectively. The superscripts ^ and ^ denote clockwise and anti-clockwise 
rotation, respectively. 

Setting the flrst variation to zero yields the following set of equations: 

"Pi,; V X B = /iiB, Pi = const., (5) 

X^; n ■ B = 0, [[Pi + 1/25']] = 0, (6) 

V; V X B = 0, V ■ B = 0, (7) 

W; n ■ B = 0, (8) 

where n is a unit vector normal to the plasma interface Xj, and [[x\] = x^+i — Xi denotes 
the change in quantity x across the interface Xj. The boundary conditions, n-B = 0, arise 
because each interface and the conducting wall is assumed to have inflnite conductivity. 
In turn, these imply the toroidal flux in each plasma region (and the poloidal flux in 
the vacuum) is constant during relaxation. Given the vessel with boundary W, the 
interfaces Xj, and the magnetic held B, Eqs. (5)-(8) constitute a boundary problem for 
the plasma pressure Pi in each region TZi. 

Minimizing the second variation subject to the constraint of the positive deflnite 
normalization = fj d'^a\C,i\'^ yields the following set of equations for the variation 
in the magnetic field b = 6B: 

Vi;Wxh = fiih, (9) 

X, ; C[[B ■ b]] + m[B{n ■ V)P]] - A^^. = 0, (10) 

n • b,,i+i = Bi^i+i ■ + ^iii ■ V X (n X B^^^+i), (11) 

V;Vxb = 0, V-b = 0, (12) 

W;n-b = 0. (13) 



Here is the normal displacement of the interface It, and A the Lagrange multiplier 
of the stability treatment, such that A < indicates a lower energy state is available. 
Using Eqs. (9)-(13) the perturbed flux through each region can be found. With a 
suitable Fourier decomposition chosen, Eq. (11) solves for the unknown coefficients of 
the perturbed field in each region. With substitution, Eq. (10) then becomes a linear 
eigenvalue equation for A. 



2.2. Tearing mode treatment 

A starting point for the treatment of tearing modes is the set of MHD equations: 

0, (14) 

J X B - Vp, (15) 
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being the fluid equation of motion, mass continuity, the adiabatic equation of state. 
Ohm's law, Faraday's law, Ampere's law, and the magnetic mono-pole condition, 
respectively. The plasma parameters change across each interface, and across surfaces 
resonant with perturbations of a given helicity. 

We solve for the plasma parameters for a zero flow plasma {i.e. v = 0) in "outer" 
regions away from the resonant surfaces where the effects of resistivity are negligible. To 
solve, the field is written B = V^x h+gfh, where g and x are scalar functions of position 
and time and h is the helical wave-field vector. Next, x ^^nd g are expanded as a Fourier 
perturbation, and solutions to the linearised Beltrami equation found. The ODE for 
Xi{f), the radial envelope of the linear Fourier perturbation for x, integrates to a jump 
condition in Xi('^)Vxi('^) at each interface, expressed in terms of equilibrium parameters. 
The plasma growth rate, obtained by linearising Faraday's law and substituting for 
E as determined by Ohm's law, is proportional to A' = [xi(r)'/xi(r)]|^l, such that 
A' = denotes marginal stability, and A' > instability. The final expression for A' 
is a function of the equilibrium parameters in the resonant region, as well as jumps in 
equilibrium parameters across the interfaces. 



3. MRXMHD cylindrical plasmas 

Solutions in an azimuthally-symmetric, axially-periodic cylinder (with axial periodicity 
length L = 2'kR) are available in Hole et al [1]. In the cylindrical co-ordinate system 



(r, 6, z) they are: 

Vi : B = {0, kiJi{fiir), kiJo{fiir) }, 

Vi : B = {0, kiJi{fiir) + diYi{fiir), kiJo{fiir) + diVo^fiir) }, (21) 

V :B = {0, BX/r, B^ }, 

where ki, di G 9?, and Jq, -'^i and Yq, Yi are Bessel functions of the first kind of order 0, 
1, and second kind of order 0, 1, respectively. The terms B'q and -Bj are constants. 
The constant di is zero in the plasma core Pi, because the Bessel functions Yq^hit) 
and Yi(/iir) have a simple pole at r = 0. Radius is normalized to the plasma- vacuum 
boundary, located at r = 1. The equilibrium is constrained by the 4iV + 1 parameters: 

{/ci, kN, d2, dN, fii, /Uat, ri, ...rAr_i, r^, Bg, B^}, (22) 

where are the radial positions of the ideal MHD barriers, and is the radial 
position of the conducting wall. Equivalent representations, and the mapping between 
these solutions has been discussed in earlier work [1-3]. 

3.1. Stability from a variational principle 

We have assessed stability using a Fourier decomposition in the poloidal and axial 
directions for the perturbed field b = V x a and the displacements of each interface. 
That is, 

b = be*('"^+^") , = Xie^('"^+'^^) , (23) 

where m, k are the Fourier poloidal mode-number and axial wave-number, and b and 
Xi are complex Fourier amplitudes. Under these substitutions, and after solving for the 
field in each region, Eq. (10) reduces to an eigenvalue matrix equation 77 ■ X = AX 
with column eigenvector X = (,^1, , eigenvalue A, and rj a N x N tridiagonal real 

matrix. 

3.2. Tearing mode stability 

In the helical coordinate u = mO + kz, a divergence-less B can be written 

B(r, u) = Vx('^, ^i) X h + 5((r, u)h, (24) 

where x is a the helical flux, and g a helical field. The vector h is defined by 
h = /(r)Vr X Vm, where /(r) = r/{m? + A;^r^) is a metric term. As in Tassi et al 
we search for helical perturbations of the form 

X(r, n, t) = xo(r) + Xi(r)e^*+^", g{r, n, t) = go{r) + gi{r)e^'^'\ (25) 

In this representation, resonant surfaces are those for which Xo(^) = 0- The equilibrium 

field satisfies the Beltrami equation, giving rise to = fi'o(^)/Xo(^)i such that the 

rotational transform can be written 

R rK^o(r)/xo(r) + m 

h = X . (26) 

r mgo{r)Xo{r)-rK 



By writing the incompressible velocity field in a similar form to Eq. (24), and expanding 
continuity to first order, it is possible to show perturbations in the fiow, pressure and 
mass density do not aff'ect marginal stability. 

In each of the plasma regions, projections of the linearised Beltrami equation along 
h and Vr yield 

9i = g'^{r)/^^{r)xi{r), (27) 
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where /i' vanishes everywhere except at r = Vstep-, where it becomes singular. These 
are identical to Eqs. (27) and (29) of Tassi et al . Equation (28) reduces to a 
Bessel or modified Bessel equation in the transformed variables xi = — K^|r 

and xy = \K\r. In the z'th region, and either side of the resonant surface, solutions 
are different combinations of Bessel or modified Bessel functions with undetermined 
coefficients Ci,Ai, Q^sy-^i-s and Ci+s^Aj+s, respectively. As only the ratio Xi{r)/xi{r) 
appears in A', its value is unaffected by normalizing C = 1 in each interval and region. 
The requirement of boundedness at r = 0, and the presence of perfectly conducting wall 
implies 

Xi(0) = Xo, XiK) = 0. (29) 

Noting that the perturbed fiux must be continuous across each interface, Eq. (28) can 
then be integrated about each interface to yield 
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(26) to eliminate ^rfX. Solving equilibrium for Xoif) gives 



The parametric dependence can also be examined by solving for Xo{r) and using Eq. 

9o{r) 
X'oir) ■ 

Xoir) = B 
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Finally, eliminating Xo('^) and -^tt^, Eq. (30) can then be rewritten 
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where 



G{B, /i, r, R, m, k) = rfiF{B, r, R, m, k) 



mR/r — KRhr/R 



(32) 
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hm + kR 

With C = 1 everywhere, the tearing mode parameter A' is a function of A in each 
interval, which is uniquely determined by the set of constraints given by Eq. (32) at 
each interface. That is, the inner boundary condition (28) yields Ai = 0. If for instance 
the resonant surface lies in 7?,2, Eq. (32) evaluated at interface Xi and X2 solves for 
in terms of Ai and h^2+s in terms of Ay, respectively. The conducting wall boundary 
condition solves for Ay. 



Changes in the field strength B at any interface enter Eq. (30) through the solution 
to Xoi'r), given by Eq. (31). Stability is hence a property of the rotational transform, 
the position of the barriers, the Lagrange multipliers, the rotational transform, and any 
jumps in the pressure or rotational transform across the interfaces. Our working reduces 
to Tassi et al in the limit of no pressure, field or rotational transform jumps across the 
interfaces, and no vacuum. 

4. Stability for an RFP-like configuration 

We have compared stability conclusions using variational and tearing mode treatments 
for an illustrative configuration. The example chosen is guided by earlier detailed 
working [1], where the core Lagrange multiplier was /xi = 2. The first interface is 
placed at r = 0.5, and the axial periodicity length chosen to be 207r, such that the 
effective aspect ratio r/R = 1/10 is small. The jump in safety factor between the 
internal interface and the plasma-vacuum boundary has been chosen to resemble Hole 
et al , subject to the different R values used for the two treatments {R = l/(27r) in Hole 
et al). We have used = 3.6, which requires ^2/^2 = 0.77 for the rotational transform 
profile to be continuous. A second motivation for this choice is the similarity to q profiles 
of high confinement reverse field pinches, such as the Madison Symmetric Torus [12] and 
RFX-mod [13], although the change in n is greatly exaggerated in this work. The plasma 
pressure is selected by the parametrization (3i = pi/{By /2fiQ), P2 = P2/{By/2fio). 
Except for the final scan over a pressure-less plasma is assumed {i.e. A = /32 = 0). 
Figure 1(b) shows the q profile for the chosen equilibrium, where q = 1/ i. 

Figure 2(a) is a dispersion curve for m = 1 modes, showing A computed using 
the variational treatment, and — r^A' computed for modes resonant within the plasma. 
Marginal stability corresponds to A = and A' = 0: these overlap identically in Fig. 2 (a). 
Modes with n = —16 and n = 12 corresponds to a perturbation near-resonant with the 
outer and inner interfaces (g(r2) ~ —1/16 and g(ri) ^ 1/12) respectively. 

In our variational treatment, we have prescribed no relationship between ^ and b 
in the relaxed regions. As such, excepting at the ideal interfaces, field line resonance in 
such plasmas is not explicitly resolved. Expressions can however be constructed which 
provide an estimate of the localization of the mode r^ efr, and a convenient choice is 
'^seff = J2iLi{''^i-^iy J where the eigenvectors are normalised such that X ■ = 1, with 
H the Hermitian. 

Figure 2(b) shows a comparison of off to r^, in which modes unstable to variational 
and tearing modes have been identified. Agreement between es and is qualitatively 
good in the interval over which the plasma in unstable, and excellent near the interfaces. 
The n = —16 and n = 9 modes are near resonant with the outer and inner interface, 
respectively. As shown in Fig. 3(a), a stability scan with wall radius indicates that 
in the limit — 1, modes for n < are wall-stabilized. In the limit that the outer 
interface is made resonant with the n = — 16 tearing mode (for example by changing 
i?). A' ^ 00. This mode is the current driven external kink of ideal MHD. Conversely, 
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Figure 2. Dispersion curves (a) and mode localization (b) of m = 1 modes of a 
pressure-less MRXMHD plasma with q profile given by Fig. 1(b). In panel (a), 
the solid lines are eigenvalues (A) of the MRXMHD treatment, and represent 
different eigenfunctions, while the points are values of —0.002 x A' determined 
from tearing mode analysis of Sec. 3.2. The vertical dashed lines correspond to 
zeros in A. Panel (b) shows the resonant surfaces of tearing modes (points), 
and effective localization ?^s,e// of modes using the variational treatment (solid 
line). The solid points and cross-hairs denote stable and unstable tearing modes, 
respectively. The heavy solid line denotes solutions for which A < 0, and the 
dashed vertical lines correspond to marginal stability, A = 0. 



the unstable range for n > 12 is only very weakly affected by the wall position. If the 
inner interface is mode resonant with the perturbation, A' ^ oo, and the mode is ideal 
unstable. This is the internal kink of ideal MHD. 

Recently, Mills et al [6] demonstrated that one can unify ideal and relaxed 
variational treatments by extending the relationship between b and ^ through the 
Newcomb gauge a = ^ x B. If this variational treatment is followed, rational surfaces 
do explicitly enter the expression for ^ as derived from br, the variation in the radial 
part of the magnetic field, and so both relaxed and tearing modes become localised at 
the resonant surface r^. 

The findings of Fig. 2, obtained for a pressure-less plasma, agree with that of 
Furth et al [14], who showed that for cylindrical pressure-less plasma with no vacuum, 
A' is linear with the second variation in the magnetic energy driving the tearing mode. 
We have also compared stability conclusions drawn from variational and tearing mode 
treatments as a function of (3. We find that while both variational and tearing mode 
treatments have the same stability trends with /3, the marginal stability limit of the 
variational treatment (for a given m and n) is lower than that of tearing modes. A 
study is ongoing into the cause of this discrepancy, as well as formally relating to 
S' 

Finally, Fig. 3(b) is a plot of the marginal stability boundary (A = 0) in n — /9 
space for m = 1 eigenmodes of the variational treatment. The two pressure profile 




Figure 3. Wall stabilization (a) and marginal stability pressure dependence 
(b) of MRXMHD plasmas unstable to m = 1 modes. Panel (a) shows the 
dispersion curves of MRXMHD plasmas unstable to m = 1 modes as a function 
of n for different conducting wall radii. In panel (b) marginal stability n — (3 
space is shown for different pressure profile configurations. The solid lines are 
for the pressure in the outer region set to zero {j32 = 0), while the dashed line 
corresponds to zero core pressure {j3i = 0). 

configurations that have been studied are /3i = and = 0. Trends in the marginal 
stability boundary can be understood by relating the radial location of the mode 
resonant surface to the analog of radial pressure gradient in the MRXMHD model: 
the sign and magnitude of nearby pressure jumps. For n > modes resonant near the 
first interface, an increasing core pressure increases the pressure drop across the first 
interface, and so destabilises the plasma. Conversely, increasing the edge pressure leads 
to a pressure jump across the first interface, and so stabilises the internal modes. For 
[32 > 0.3 all m = 1 internal modes {n > 0) are completely stabilised. For the n < 
modes resonant near the edge, changes in the core pressure have little effect, while 
increasing the edge pressure destabilises the plasma. 

5. Conclusions 

We have computed the stability of multiple relaxation region MHD (MRXMHD) plasmas 
using both a variational and a tearing mode treatment, evaluated in a periodic cylindrical 
configuration. The marginal stability conclusions of the two treatments for a zero P 
plasma, as well as the trends with /3, appear to be identical, in agreement with earlier 
analytic working by Furth et al [14]. Some discrepancy exists between the marginal 
stability boundaries of variational and tearing mode treatments for nonzero /3, with 
the stability limit of variational plasmas lower than that of tearing modes. A study is 
underway to resolve this discrepancy. 

The overarching aim of this work is to elucidate the nature of perturbations 



available to MRXMHD equilibria, which in turn, are motivated by our quest for 
mathematically rigorous solutions of ideal MHD force balance in 3D geometry. Our 
working builds of Tassi et al to nonzero (3 multiple region relaxed plasmas with vacuum, 
and complements Mills et al , who demonstrated that that one can unify ideal and 
relaxed variational treatments through the Newcomb gauge. Combined, these results 
suggest that MRXMHD configurations which are stable to ideal perturbations plus 
tearing modes are automatically in a stable equilibrium state. 

In ongoing work we are developing a faster numerical algorithm for the construction 
of 3D MRXMHD plasmas, and exploring the maximum pressure jump an interface can 
support before it is destroyed by chaos. 
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